// Copyright (C) 2011  Andrew H. Chan
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

#include <iostream>
#include <iomanip>
#include <vector>
#include <cstdlib>
#include <math.h>
#include <assert.h>
#include <algorithm>

#include "differenceCoverSample.hpp"

#include "ntSeq.hpp"
#include "LarssonSadakane.hpp"
#include "findMedian.hpp"

typedef unsigned long long uint64;

using namespace std;

// Luk, Wong: Two new quorum based algorithms for distributed mutual exclusion. 1997.
const uint64 DCS::dc0[] = {0};
const uint64 DCS::dc1[] = {0,1};
const uint64 DCS::dc2[] = {0,1,2};
const uint64 DCS::dc3[] = {0,1,2,4};
const uint64 DCS::dc4[] = {0,1,2,5,8};
const uint64 DCS::dc5[] = {0,1,2,3,7,11,19};
const uint64 *DCS::dcAll[] = {dc0,dc1,dc2,dc3,dc4,dc5};
const uint64 DCS::dcSize[] = {1,2,3,4,5,7};

DCS::DCS(ntSeq const &T, uint64 logV) : n(T.size()), logV(logV) {
    assert(logV>0);
    assert(n>0);

    v = 1<<logV;
    vMask = v-1; // compute vMask

    computeDifferenceCover(logV);
    assert(!verifyDC());

    computeDelta(); // sets d
    assert(!verifyDelta());

    computeRank();
    computeISA(T);
}

// compute difference cover
void DCS::computeDifferenceCover(uint64 logV) {
    // handle small cases
    if(logV < 6) {
        DC.resize(dcSize[logV]);
        for(uint64 ii=0; ii<dcSize[logV]; ++ii) DC[ii] = dcAll[logV][ii];
        return;
    } else { // handle the rest
        // compute r
        double tmpR = (2.*sqrt((double)(81 - 6*(13-v))) - 18.)/(24.);
        uint64 r = tmpR>0 ? (int)ceil(tmpR) : 0;
        assert(24*r*r + 36*r + 13 >= v);
        DC.resize(6*r+4);
        DC[0]=0;
        for(uint64 ii=1; ii<r+1; ++ii) DC[ii] = 1;
        DC[r+1] = r+1;
        for(uint64 ii=r+2; ii<2*r+2; ++ii) DC[ii] = 2*r+1;
        for(uint64 ii=2*r+2; ii<4*r+3; ++ii) DC[ii] = 4*r+3;
        for(uint64 ii=4*r+3; ii<5*r+4; ++ii) DC[ii] = 2*r+2;
        for(uint64 ii=5*r+4; ii<6*r+4; ++ii) DC[ii] = 1;
        for(uint64 ii=1; ii<6*r+4;++ii) DC[ii] = DC[ii-1]+DC[ii];
        return;
    }
}

void DCS::computeDelta() {
    d.resize(v);
    for(uint64 ii=0; ii<DC.size(); ++ii)
        for(uint64 jj=0; jj<DC.size(); ++jj)
            d[modDiff(DC[jj],DC[ii])] = DC[ii];
}

void DCS::computeRank() {
    rank.resize(v);
    rank[0] = 0;
    uint64 eltRank = 0;
    for(uint64 ii=1; ii<v; ++ii) {
        if(eltRank < DC.size()-1 && ii==DC[eltRank+1]) ++eltRank;
        rank[ii] = eltRank;
    }
    h = rank[n%v];
}

void DCS::computeISA(ntSeq const &T) {
    // compute Dn
    uint64 lenDn = DC.size() * (n/v) + h + 1;
    vector<uint64> Dn;
    Dn.reserve(lenDn);
    for(uint64 ii=0; ii<DC.size(); ++ii) {
        uint64 maxjj = n/v + (DC[ii]<=n%v ? 1:0);
        for(uint64 jj=0; jj<maxjj; ++jj) Dn.push_back(DC[ii] + jj*v);
    }
    assert(Dn.size() == lenDn);

    l.resize(Dn.size());
    multiQsort(T, Dn, 0, 0, Dn.size()); // v-orders the suffixes; modifies l
    vector<uint64> ss = l; // debug
    LarssonSadakane::computeInverseSuffixArray(l);
    assert(!LarssonSadakane::verifyInverseSuffixArray(ss, l, 0));
    assert(!verifyDCS(T,ss));
}

// v-order the suffixes in Dn
// modifies l
void DCS::multiQsort(ntSeq const &T, vector<uint64> &Dn, uint64 charIndex, uint64 index, uint64 length) {
    assert(length>0); // multiQsort should not be called with length 0 
    //if(length==0) return;

    uint64 charIncr = 21*charIndex;
    // selection sort serves as the base case
    if(length < 7) { // selection sort if partition has fewer than 7 elements 
        uint64 pLow=index, pLast=index+length-1;
        while(pLow < pLast) {
            uint64 minKey = packedKey(T,Dn[pLow]+charIncr);
            uint64 pEqual = pLow+1;
            for(uint64 ii=pLow+1; ii<=pLast; ++ii) { // search for minimum element
                uint64 tmpKey = packedKey(T,Dn[ii]+charIncr);
                if(tmpKey < minKey) {
                    minKey = tmpKey;
                    swap(Dn[ii], Dn[pLow]); // swap to beginning of selection group 
                    pEqual=pLow+1; // new min found, so reset pEqual to point to one past end of min elements
                } else if(tmpKey == minKey) {
                    swap(Dn[ii], Dn[pEqual]); // swap to min elements group at beginning of selection group
                    ++pEqual;
                }
            }
            uint64 numEqual = pEqual-pLow;
            
            if(numEqual>1 && charIndex<(v-1)/21) {
                multiQsort(T, Dn, charIndex+1, pLow, numEqual); // recursively sort equal group
            } else {
                for(uint64 ii=pLow; ii<pLow+numEqual; ++ii)
                    l[mu(Dn[ii])] = pLow+1;
            }
            pLow = pEqual; // prepare to sort next selection group
        }
        if(pLow==pLast) {
            l[mu(Dn[pLast])] = pLast+1;
        }
        return;
    } else { // quicksort if partition has 7 or greater elements
        // choose pivot
        uint64 pivot;
        if(length <= 50) { // if partition has 50 or fewer elements, select the median of first, middle and last elements
            uint64 a = packedKey(T,Dn[index]+charIncr), b = packedKey(T,Dn[index+(length>>1)]+charIncr), c = packedKey(T,Dn[index+length-1]+charIncr);
            pivot = a > b ? (b > c ? b : (a > c ? c : a)) : (a > c ? a : (b > c ? c : b));
        } else { // if partition has more than 50 elements, select median of 5 evenly spaced elements
            uint64 incr = length>>2;
            pivot = findMedian5(packedKey(T,Dn[index]), packedKey(T,Dn[index+incr]+charIncr), packedKey(T,Dn[index+incr+incr]+charIncr), packedKey(T,Dn[index+incr+incr+incr]+charIncr), packedKey(T,Dn[index+length-1]+charIncr));
        }
        uint64 pCurLow = index, pLow = index, pCurHigh = index+length-1, pHigh = index+length-1;
        uint64 curKey;
        while(true) {
            while(pCurLow <= pCurHigh && (curKey=packedKey(T,Dn[pCurLow]+charIncr)) <= pivot) {
                if(curKey == pivot) swap(Dn[pCurLow],Dn[pLow++]); // swap equal elements to beginning
                ++pCurLow;
            }
            while(pCurHigh > pCurLow && (curKey=packedKey(T,Dn[pCurHigh]+charIncr)) >= pivot) {
                if(curKey == pivot) swap(Dn[pCurHigh],Dn[pHigh--]); // swap equal elements to end
                --pCurHigh;
            }
            if(pCurLow>=pCurHigh) break;
            swap(Dn[pCurLow++],Dn[pCurHigh--]); // swap low with high
        }
        pCurHigh = pCurLow-1;
        // pCurLow points one beyond the smaller elements
        // pCurHigh points one to the left of the larger elements

        // swap equal elements from left toward center
        uint64 numEqLeft = pLow-index; // number of elements equal to pivot
        uint64 numLow = pCurLow-pLow; // number of elements less than pivot
        if(numEqLeft > numLow)
            for(uint64 ii=index, jj=pLow, numSwap=numLow; numSwap>0; --numSwap, ++ii, ++jj) swap(Dn[ii], Dn[jj]);
        else
            for(uint64 ii=index, jj=pCurLow-numEqLeft, numSwap=numEqLeft; numSwap>0; --numSwap, ++ii, ++jj) swap(Dn[ii], Dn[jj]);

        // swap equal elements from right toward center
        uint64 numEqRight = index+length-1-pHigh; // number of elements equal to pivot
        uint64 numHigh = pHigh-pCurHigh; // number of elements greater than pivot
        if(numEqRight > numHigh)
            for(uint64 ii=pCurLow, jj=index+length-numHigh, numSwap=numHigh; numSwap>0; --numSwap, ++ii, ++jj) swap(Dn[ii], Dn[jj]);
        else
            for(uint64 ii=pCurLow, jj=index+length-numEqRight, numSwap=numEqRight; numSwap>0; --numSwap, ++ii, ++jj) swap(Dn[ii], Dn[jj]);

        uint64 numEqual = numEqLeft + numEqRight;

        if(numHigh > numEqual && numHigh > numLow) { // larger is largest partition
            if(numLow>0) multiQsort(T, Dn, charIndex, index, numLow); // sort smaller
            if(numEqual>0 && charIndex<(v-1)/21) {
                multiQsort(T, Dn, charIndex+1, index+numLow, numEqual); // sort equal
            } else {
                for(uint64 ii=index+numLow; ii<index+numLow+numEqual; ++ii)
                    l[mu(Dn[ii])] = index+numLow+1;
            }
            if(numHigh>0) return multiQsort(T, Dn, charIndex, index+length-numHigh, numHigh); // sort larger
        } else if(numLow > numEqual && numLow > numHigh) { // smaller is largest partition
            if(numHigh>0) multiQsort(T, Dn, charIndex, index+length-numHigh, numHigh); // sort larger
            if(numEqual>0 && charIndex<(v-1)/21) {
                multiQsort(T, Dn, charIndex+1, index+numLow, numEqual); // sort equal
            } else {
                for(uint64 ii=index+numLow; ii<index+numLow+numEqual; ++ii)
                    l[mu(Dn[ii])] = index+numLow+1;
            }
            if(numLow>0) return multiQsort(T, Dn, charIndex, index, numLow); // sort smaller
        } else { // equal is largest partition
            if(numHigh>0) multiQsort(T, Dn, charIndex, index+length-numHigh, numHigh); // sort larger
            if(numLow>0) multiQsort(T, Dn, charIndex, index, numLow); // sort smaller
            if(numEqual>0 && charIndex<(v-1)/21) {
                return multiQsort(T, Dn, charIndex+1, index+numLow, numEqual); // sort equal
            } else {
                for(uint64 ii=index+numLow; ii<index+numLow+numEqual; ++ii)
                    l[mu(Dn[ii])] = index+numLow+1;
            }
        }
        return;
    }
    assert(false);
}

// verify difference cover
int DCS::verifyDC() const {
    cout << "Verifying difference cover for v: " << v << endl;
    assert(DC.size()>0);
    vector<bool> check(v,false);
    for(uint64 ii=0; ii<DC.size(); ++ii)
        for(uint64 jj=0; jj<DC.size(); ++jj)
            check[modDiff(DC[ii],DC[jj])] = true;
    for(uint64 ii=0; ii<v; ++ii)
        if(!check[ii])
            return 1;
    return 0;
}

int DCS::verifyDelta() const {
    cout << "Verifying delta for v: " << v << endl;
    vector<bool> covered(v,false);
    for(uint64 ii=0; ii<DC.size(); ++ii) covered[DC[ii]] = true;
    for(uint64 ii=0; ii<v; ++ii) {
        for(uint64 jj=0; jj<v; ++jj) {
            uint64 k = delta(ii,jj);
            if(!covered[(ii+k)%v] || !covered[(jj+k)%v]) {
                return 1;
            }
        }
    }
    return 0;
}

int DCS::verifyDCS(ntSeq const &T, vector<uint64> const &ss) const {
    cout << "Verifying DCS for v: " << v << endl;
    for(uint64 ii=0; ii<T.size(); ++ii) {
        for(uint64 jj=ii+1; jj<T.size(); ++jj) {
            // compute lcp
            string suffixA = T.getSuffix(ii);
            string suffixB = T.getSuffix(jj);
            uint64 kk=0;
            for(; kk<v-1; ++kk) {
                if(suffixA[kk] != suffixB[kk]) break;
            }
            assert(kk<=v-1);
            if(kk==v-1) {
                if(suffixA < suffixB) {
                    if(!lt(ii,jj)) {
                        return 1;
                    }
                } else {
                    if(!lt(jj,ii)) {
                        return 1;
                    }
                }
            }
        }
    }
    return 0;
}
